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Abstract 

We consider the effects of parameter uncertainty on the optimal radiation schedule in the 
context of the linear-quadratic model. Our interest arises from the observation that if inter¬ 
patient variations in normal tissues and tumor sensitivities to radiation or sparing factor of the 
organ at risk (OAR) are not accounted for during radiation scheduling, the performance of the 
therapy may be strongly degraded or the OAR may receive a substantially larger dose than 
the maximum threshold. This paper proposes two radiation scheduling concepts to incorporate 
inter-patient variability into the scheduling optimization problem. The first approach is a robust 
formulation that formulates the problem as a conservative model that optimizes the worst case 
dose scheduling that may occur, assuming that the parameters vary within given intervals. The 
second method is a probabilistic approach, where the model parameters are given by a set of 
random variables. Our probabilistic formulation insures that our constraints are satisfied with 
a given probability, and that our objective function achieves a desired level with a stated prob¬ 
ability. We used the same transformation as [37] to reduce the resulting optimization problem 
to two dimensions. We showed that the optimal solution in the absence of uncertainty in the 
tumor radio-sensitivity parameters (a and (3) occurs at one of the corners of the feasible region. 
However if we incorporate uncertainty in a and (3 into the optimization problem, this result 
does not hold anymore. In this case, we showed that the optimal solution lies on the boundary 
of the feasible region and we implemented a branch and bound algorithm to find the global 
optimal solution. We demonstrated how the configuration of optimal schedules in the presence 
of uncertainty compares to optimal schedules in the absence of uncertainty (conventional sched¬ 
ule). We observed that if the number of fractions in the optimal conventional schedule is the 
same as the robust and stochastic solutions, it is preferable to administer equal or smaller total 
dose. In addition if there exist more (fewer) treatment sessions in the probabilistic or robust 
solution compared to the conventional schedule, a reduction in total dose squared (total dose) 
will be expected. Finally, we performed numerical experiments in the setting of head-and-neck 
tumors including several normal tissues to reveal the effect of parameter uncertainty on optimal 
schedules and to evaluate the sensitivity of the model to the choice of key model parameters. 
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1 Introduction 


The building block for virtually all mathematical models of radiation response is the linear-quadratic 
model (LQ), which matches well with experimental data across a wide range of clinically relevant 
radiation doses and fractionation schemes m and 0) The basic model states that if a collection 
of cells is exposed to N fractions of radiation, di Gy (SI derived unit of ionizing radiation) in i th 
fraction, the reproductively viable fraction of cells after the exposure is given by e~ ad i+Pd ?. 

The two parameters a and (3 depend on the specific tissue that is being irradiated. The parameter 
a represents killing of cells from a single track of radiation, and (3 represents the killing of a cell 
via two independent tracks of radiation 023 - There are several mathematical extensions to the 
LQ framework to incorporate additional biological phenomena such as repopulation of the tumor 
population between fractions, re-oxygenation of the tumor (this is required for some radiation 
therapy to be effective), the effectiveness of DNA repair mechanisms between fractions, and the 
redistribution of tumor cells within the cell cycle. Taken together these four extensions are often 
referred to as the ‘4Rs’ and there have been several works based on these extensions [451] . 

When radiotherapy is used in the clinical setting, it is necessary to ensure that the treatment 
avoids excessive toxicity in normal tissues in the vicinity of the tumor. Therefore it is necessary to 
ensure that the radiation absorbed by the surrounding normal tissue falls within desired constraints. 
Hence the ultimate goal in radiotherapy is maximizing tumor damage while ensuring that the level 
of normal tissue toxicity does not exceed a given threshold. The standard approach for measuring 
tumor damage and tissue toxicity is via the linear quadratic model and the biologically equivalent 
dose (BED), respectively ([15] and USD- 

Most radiation treatments are currently administered in equal fractions five days a week, for 
6 weeks total. Over the past few decades there have been several mathematical works that have 
studied the survival benefit of various fractionation schedules for a wide range of cancers. In [6], 
that most other radiobiological make similar time-dose predictions as the LQ formalism. In that 
work they used the LQ model in combination with Lea-Catcheside time factor, which takes into 
account dose protraction or fractionation and DNA repair between fractions. Yang and Xing m) 
explored the influence of the ‘4Rs’ of radiobiology on external beam radiotherapy for fast and 
slowly proliferating tumors and conclude that including repair effects in the BED model may give 
rise to optimal non-uniform fractionation schedules. Mizuta et al. (E3) presented a mathematical 
model that minimizes the radiation effect on the late responding normal tissues while keeping the 
effect of radiation on the tumor constant. They showed that the multi-fractionated irradiation 
with a constant dose is better if the ratio of ( % ) / (%) is less than the ratio of 

the dose received by the normal tissue, while Hypo-Fractionated irradiation is better otherwise. 
Unkelbach et al. (|43|) studied the interdependence of the optimal fractionation scheme and the 
spatial dose distribution in the normal tissues. In particular, they derived a criterion under which 
a Hypo-Fractionated regimen is indicated for both parallel and serial OARs. In a very recent work 
[36], a formulation of the optimal fractionation problem that includes multiple normal tissues has 
been considered. They established sufficient conditions under which equal-dosage or single-dosage 
fractionation is optimal. In recent work m) the authors investigated optimal fractionation for a 
mouse model of glioblastoma, in this work they found that non-standard fractionation schedules 
lead to improved survival times; a finding that was verified in experimental studies. In [2| this work 
was extended to include a richer set of toxicity constraints. 

Until very recently, there was no work that precisely described the optimal fractionation sizes in 
the presence of multiple normal tissues. In particular, most works considered the optimal schedule 
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with respect to a single normal tissue. However in practice, there are usually at least two healthy 
structures in the vicinity of the tumor. Saberian et al. considered several normal tissues in their 
study, however they were unable to find the closed form solution to the problem for all possible cases 
and they only discussed the sufficient conditions under which equal-dosage fractionation is optimal 
f36 j. In p] two simultaneous normal tissue toxicity constraints were implemented. In a very recent 
work, Saberian et al. m found the closed form solution to the problem of optimal fractionation 
while maintaining multiple simultaneous normal tissue constraints without considering any pre¬ 
sumptions about the configuration of the optimal solution. They solved the problem to optimality 
by instead solving a two-variable linear program with two additional nonlinear constraints. 

An important result emerging from recent work is that the sparing factor of normal tissues 
and the magnitude of the a/(5 ratio for both normal tissues and the tumor determine the optimal 
radiation schedule (0, m, eh and m)- Therefore the optimal fractionation schedule is acutely 
sensitive to perturbations in these parameters. One consequence of this sensitivity is the following: 
an optimal fractionation schedule will have been derived for a fixed set of parameter values (called 
the nominal values), but for a specific patient with a distinctly different set of parameter values 
this schedule is no longer optimal, and in fact may have poor performance. The uncertainties in 
radiotherapy treatment can be categorized into two groups: geometric and inter-patient variabil¬ 
ity. Target volumes take account of geometric uncertainties such as organ motion, inaccuracies 
or variations in treatment set-up, patient positioning errors and fluctuations in machine output. 
Several studies addressed these uncertainties using different techniques. Stroom et al [3Hj devel¬ 
oped a method for the automatic calculation of planning target volume margins as a means for 
incorporating geometric uncertainties in the region that is irradiated. The traditional approach to 
dealing with uncertainty in IMRT (considering a margin surrounding the tumor volume) increases 
the radiation exposure of healthy tissue. Chan et al. |8j developed a robust framework to incorpo¬ 
rate uncertainty in the probability distribution that describes breathing motion, and showed that 
a treatment plan obtained from the robust formulation delivers 38% less dose to the OARs than 
the traditional solution, while providing the same level of protection against breathing uncertainty. 
In a similar work, Chu et al. [[9] used a robust optimization approach to find the IMRT treatment 
plans while considering patient motion and setup uncertainties. More specifically they included 
uncertain voxel location in their model and as a consequence the delivered dose became a random 
variable. They designed a mathematical model constructing plans that are more adept at sparing 
healthy tissue while maintaining the prescribed dose to the target under uncertainty. In [02], two 
methods to account for range uncertainties, one method using a probabilistic approach and the 
other applying methods from robust linear programming were presented to find optimized treat¬ 
ment plans for intensity modulated proton therapy. Both methods greatly reduced the sensitivity 
to range uncertainties of the resulting treatment plans. A modification of the worst case optimiza¬ 
tion was applied to a clinical case by Pflugfelder et al. [32]. In addition to the robust optimization, 
stochastic programming has also been used to account for organ motion and setup errors in IMRT 
optimization (see [22], [23]), e.g. Unkelbach developed a planning method that accounts for the 
probabilistic dwelltime of a tumor evaluated from multiple CT scans [43J. 

Inter-patient variability is due to heterogeneity in patient-specific variables such as the sen¬ 
sitivity of their normal tissues and tumor to radiation, and the growth rate of their tumor. In 
several cancers there have been multiple subtypes discovered driven by distinct genetic pathways 
and having distinct phenotypic behaviors such as growth parameters and response to therapy (e.g. 
glioblastoma [28], breast cancer [29], head and neck cancer [12], melanoma [26] and many others). 
A distinct possibility is that there is still significant patient variability within these subtypes. In 
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fact this inter patient heterogeneity is a large reason for the pursuit of personalized medicine [18]. 
Given current technologies it is difficult to measure tumor response parameters a and ft during 
treatment due to confounding effects such as protracted cell death (T3j, cell cycle arrest |4], and 
radiotherapy mediated immune response m • Furthermore, toxicity effects often do not show up 
until several months or even years after conclusion of therapy, and it is therefore not possible to 
learn the tumor response properties of normal tissues during treatment. In this paper we concen¬ 
trate on the modeling the uncertainties arising in inter-patient variations, which we will use later 
in an optimization method for radiotherapy scheduling. In a recent unpublished manuscript [I] 
(appearing on web after the first version of current manuscript), the authors developed a similar 
model where the uncertain parameters is assumed to take values in a given interval. Unfortunately 
this optimization method may lead to an overly pessimistic solution, and furthermore they did not 
include the uncertainty in the tumor radio-sensitivity parameters (a and j3). To the best of our 
knowledge, the present study is the first to address the stochastic and linear uncertainty generated 
by inter patient heterogeneity. 

We present a mathematical formulation of the optimal fractionation problem in the presence of 
multiple normal tissues incorporating uncertainties in model parameters based on the LQ model 
adjusted for tumor proliferation with a time lag. This formulation allows for the parametric uncer¬ 
tainty to take two forms. First a minimal underlying stochastic model of the uncertain parameters 
is assumed to be known and every parameter, independently of other entries, takes values in a given 
interval. We formulate our problem as a model whose solution must be feasible for all realizations 
of the parameters, and even a small violation of the constraints cannot be tolerated. This method 
may lead to an overly pessimistic solution, therefore we develop additional models where we assume 
that the uncertain parameters are characterized by a probability distribution and we reformulate 
our optimization problem to now insure that our constraints are satisfied with a given probability, 
and that our objective function achieves a desired level with a given probability. We examine the 
mathematical properties of the optimal fractionation scheme in various models. The results are 
discussed in the context of head and neck tumors. As a generalization, we broadly consider the 
effects of parametric uncertainty on the structure of optimal fractionation schedules. 

The organization of the remainder of the paper is as follows. In section two we describe the 
problem formulation in the setting of fixed parameters (the nominal setting), and then formulate 
the robust and probabilistic counter parts of this nominal problem for various uncertainty sets. In 
the next section we describe our solution methods for the problems presented in section two. In 
section four we solve our optimization problems for the specific case of head and neck carcinomas 
and generalize the results to the other tumors. In section five we summarize our results and discuss 
the implications of our findings. 

2 Model of uncertainty and robust formulation 

In this section, we first define an objective function derived from the standard linear-quadratic 
model of radiotherapy response. We next discuss the constraints that are present in our optimization 
problem, which are derived by maintaining a fixed level of normal tissue damage for a variety of 
tissue types. Finally we incorporate parameter uncertainty by formulating robust and probabilistic 
versions of our optimization problem. 
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2.1 The nominal formulation 


We now consider the problem of finding fractionation schedules that lead to maximal tumor re¬ 
duction while maintaining acceptable levels of normal tissue damage. The basic linear quadratic 
(LQ) model states that if a collection of tumor cells are exposed to N fractions of radiation with 
dj Gy (SI derived unit of ionizing radiation) in j th fraction, the reproductively viable fraction of 

cells is given by e j p j . However, reproductively viable tumor cells will eventually begin 

to reproduce, and thus the total surviving fraction of surviving cells is often adjusted to take into 
consideration the reproduction of tumor cells. A common way to model the repopulation effect 
is to assume an exponential repopulation process, see e.g., [30]. Thus the net surviving fraction 
( S) due to combined effects of radiation and repopulation after the conclusion of a fractionated 
radiotherapy treatment is given by 


S = 


— P - Efel adj+pd 2 - 


ln(2 )(T r -T k ) + 




where T r , T e and T are respectively radiation delivery duration, effective cellular doubling time 
and kick-off time (or lag before exponential growth begins). The expression ( T r — T),) + is defined 
as max(0, T r — T^). Throughout this paper, we made two important assumptions. First in order to 
consider the impact of working hour constraints on the objective function, we assume that working 
hour constraints require that radiation can only be delivered hourly between 8 am and 8 pm and 
five days per week. Second, we assume for every schedule there exist n daily fractions with equal 
time elapsed between consecutive fractions. If we define a and r as the quotient and remainder of 
^, respectively, and a' and r' as the quotient and remainder of |, respectively, we can compute T r 

8 I 12 r ~ ^ 

when a/0. When a = 0, we simply have T r = — 24 "~ 1 


as ( 2.1 


T r = 


7 o' + r', 

7 (a 1 - 1) + 5, 

7 a' + r' + ^p, 
7(a' - 1) + 5 + 


r = 0 , r' 7 ^ 0 
r = 0 , r' = 0 

r > iy / 0 

r > 1 , r 1 = 0 


( 2 . 1 ) 


A natural risk associated with radiotherapy is damage to normal tissue near the tumor. A 
further complication to this toxicity is that in any radiotherapy treatment there are often a large 
number of normal tissues exposed to radiation. In addition to the existence of a large number of 
normal parenchymal cells in the clinical target volume of the respective organ, all tumor volume 
contains various stromal tissues (e.g. blood vessels and normal connective tissue). In all these 
normal cells and structures, radiation side-effects may be different (e.g. see the effect on radiation 
on parallel, serial and dose-volume organs in |19|). A common measure of toxicity for various 
normal tissues is the biologically equivalent dose or BED HU- In particular, assume that for a 
specific normal tissue of interest the radiation response is characterized by parameters a* and /%, 
furthermore assume that this tissue is exposed to N fractions of sizes {di,..., d at} respectively, and 
lastly assume that for fraction j normal tissue is only exposed to didj Gy of radiation for a sparing 
factor 5i G (0,1]. For each normal tissue we define the maximal toxicity 


BEDT 


Di + 5i 


frDl 

OLi Ni 


where Di and Ni are tissue specific parameters. Note that Di is a tissue specific parameter, e.g. it 
is defined as maximum total dose for the serial normal tissues and maximum mean dose in parallel 
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normal tissues [36]. The toxicity constraints for all the OAR is then given by 


N 


J2(dj + —did)) < BED™ ax for 1 <i< M, 

OLi 

3 =1 


( 2 . 2 ) 


where M is the number of different normal tissues under consideration. By taking the natural 


logarithm of objective function and using (2.2) to model acceptable normal tissue damage, the 


nominal problem of finding fractionation schedules that lead to maximum tumor damage while 
maintaining acceptable levels of normal tissue damage can be modeled as 


N 


max adn + pdR — q(N) 

>n J J ' ' 


d,->0,iVeZ+ . , 

j=l 


(2.3) 


subject to 


N 


+ -did]) < Di + di-^d-, i = 1,..., M 

3 =1 


where g(N) = M 2 )^ T=] + - 


2.2 Modeling uncertainty in radiobiologic parameters 


In order to solve the optimization problem ( |2.3[ ) it is vital to know the parameters a, P, Pi/on, and di 
since the optimal fractionation schedule will depend on their value ( 123,123 and 03]). However, it is 
quite difficult to obtain accurate measurements of these parameters in a clinical setting and precise 
estimates of these values are very difficult to find. Furthermore, due to inter-patient heterogeneity 
it is possible that a wide range of parameter values are possible across the patient population. For 
example in several cancers there are a multitude of possible mutational pathways responsible for 
the creation of the tumor, e.g., breast, glioblastoma, and head & neck. As a result of this situation 


we are investigate the effect of parametric uncertainty on the solution to problem (2.3). 


We assume uncertainties presented in LQ model can take two forms: (i) estimation errors for 
parameters of constant but unknown value, and (ii) stochasticity of random variables. In the first 
case only the range of the uncertain parameters is known, specifically, we assume parameter a 
belongs to a symmetric interval [a — l a ,a + l a \ centered at a and for the second scenario we consider 
a as a continuous random variable with probability density function /. In the second case, we are 
interested in finding the optimized radiotherapy delivery schedule based on two principles: first the 
nominal values of sensitive parameters are inaccurate and we only know that they lie in a given 
and second, using the range alone may lead to an excessively high level of conservativeness and the 
the objective function may suffer as a result. 


2.2.1 Non-probabilistic robust formulation 

As mentioned above the parameters a, P, {di}f£ :1 and {Pi/ai}f£ 1 are subject to uncertainty and 
may vary amongst patients. For example the values 0.33Gy and O.lOGy are frequently assumed for 
the ratio P/a for late responding normal tissue and tumor tissue respectively. However these values 
should be considered as a rough estimate as there is little evidence |19| to show that these values can 
be generalized across a wide range of human normal-tissue endpoints and tumor histologies. For the 
sparing factor d, there has been a significant amount of effort dedicated to improving the accuracy 
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and precision of radiation therapy delivery in the past decades. However there still exist sources 
of uncertainty (e.g., patient motion, organ deformation, positioning uncertainty) which make it 
impossible to achieve full precision in estimating parameters associated with organ movements in 
radiotherapy, and thus the exact value of the sparing factor 6 is often not known. 


Here our aim is to construct a robust formulation to (2.3) that is immune to realizations of 


the uncertain parameters so long as they lie within certain sets. This approach may be the only 
reasonable alternative when the parameter uncertainty is uniformly distributed, or if no distribu¬ 
tional information is available. First we consider the case that the values of parameters a and /3 
are known and we only know that the parameters {/3i/a-i}f£i and lie in given intervals. 

These uncertainty sources can have a detrimental effect on configurations or feasibility of optimal 
schedules. In this paper, we assume a fixed amount of radiation are delivered to tumor, however 
the fraction of radiation absorbed by normal tissues is subject to uncertainty. Here we assume for 
the i th normal tissue, Pi/ai and 5i are modeled as symmetric and bounded random variables that 
take lie in given intervals, or equivalently 


— 6 [-4,-b h], and <5, E [dj — lg f , 5i + lg J, i £ {1,..., M}. 

(Xi Qti Oi% 


(2.4) 


Formally, the robust counterpart of (2.3) considering uncertainties defined in (2.4) can be written 
as 


subject to 


max 

dj>o,Nez+ 


N 

adj + f3d 2 — g(N ) 

3 =1 


£>i n £>2 n • • • n B M -i n B M 


(2.5) 


where the definition of B i,..., Bm and the derivation of the robust counterpart can be found in 
the appendix. 


2.2.2 Probabilistic optimization models 


Although (2.5) provides the strongest protection against excessive toxicity in OAR, it is also the 
most conservative solution and results in less tumor cell kill than achieved by optimizing the 
nominal formulation. To address this excessive conservativeness, we control the level of flexibility 
between robustness and performance of the optimal schedule by using a probabilistic formulation 
that provides a notion of a budget of uncertainty. We view a and f3 as continuous random variables 
with joint probability density function /(•,•) and we assume that the cdf of — and <5* in the i th 


normal tissue are Fj and Gi respectively. In addition we assume that for each i ^ and <5, are 
independent of all other random variables in the model. We will require that BED in the i th 
normal tissue does not exceed some level with a high probability. This desire can be naturally 
expressed by requiring that the BED in the i th normal tissue exceeds the maximum allowable 
BED, BED™ ax , with probability at most 1 — pi, where pt is some constant close to 1, e.g., 0.95. 
Furthermore we require that optimized schedules obtained by our robust formulations result in an 
objective value which exceed level z, with probability more than or equal to p z . Computing the 
above probabilities, we can derive 


max z — g(N) 
dj>0,NeZ+,z 


( 2 . 6 ) 
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2-«E 


3 = 1 0 


subject to 

[ E ^ =1 dj [ E ^ =1 ^ f(a, (3)df3da < (1 - p z )P(a > 0, /3 > 0) 

Jo Jo 

s 1 ns 2 n---ns M -inS M - 

Derivation and the definition of S \,... ,Sm are shown in the appendix. 


3 Solution approach 

We now turn our attention to the solution of the optimization problems presented in the previous 
section. For any fixed N, the feasible regions of models described in Section 2 are compact sets 
and the objective functions are continuous. Therefore by the extreme value theorem of Weierstrass 
optima exist. 

First we use the variable transformation described in m to simplify the formulations posed in 


(2.3), (2.5), and (2.6). The simplification is to introduce the variables 


N 


N 




(3.7) 


3 = 1 


3 = 1 


We can use the results of m and conclude that there is a feasible solution for ( |3.7[ ) if and only 
if X 2 < NY and X 2 > Y. As a consequence, by adding these constraints, we guarantee that the 
optimal doses d* can be retrieved based on the solution of the adjusted formulations. 

The proliferation term, g(N), does not depend on dose. By including tumor proliferation, we 
can also optimize the number of treatment sessions N. This can be done by first introducing the 
maximum number of radiation fractions we are willing to administer in the course of radiation 
therapy, N max . The for each 1 < N < N max we solve the simplified versions of problems posed 


in (2.3) and (2.5) using transformation introduced in (3.7) and the two additional constraints 
X 2 < NY and X 2 > Y. Finally, we choose the N that maximizes the optimal biological effect on 
the tumor and return the optimal ( X *, Y *) associated with the optimal N as the global optimal 
solution. 


The feasible region defined by each B t or S t using the transformations in (3.7) can be described 
in the following simplified form 


D 2 

Y > —?-,X- 

~ Ni' 


'NY < Di 


Pi 

i N i 


D 2 D 2 

Y < 1 X + r'Y < D + r' - 

S TV,-’ + 1 S 1 A, 


i = 1,..., M 


(3.8) 


where Ci and cj are the coefficients from B l or Si. When a and /3 are fixed, for every N , we can replace 
the conditions X 2 < NY and X 2 > Y with Y > O^X and Y < 9iX (see ^7]), respectively, where for 
1 < k < N maxi di = and the pair (X^,Y^) are the coordinates of the intersection kY = X 2 with 
polygon defined by the OAR constraints given by ( |3.8[ ). Therefore in order to solve the optimization problems 
(2.5) and (2.6) with a fixed value N € {1,..., N max }, we need to specify the corners of the convex hull defined 
by the inequalities in (3.81, Y > d^X and Y < 9\X. In the appendix, Algorithm [l] describes how to construct 
the corners of the polygon sorted in increasing order of their 2-coordinates CHn = {(Xi, Yi),..., (X*,, W)}. 


The optimal solution to (2.3) and (2.5) for every N occurs at one of these corners. 


For (2.6) a more specialized approach is required. The new formulation of (2.6) is 


max z — q(N) 

X,Y>0,z 


(3.9) 



















subject to 


/(a, (3)d(3da < (1 — p z )P(a > 0, /3 > 0) 

Si n s 2 n • • • n s M -i n s M 
x 2 < N max Y 
X 2 > Y 


(3.10) 

(3.11) 

(3.12) 

(3.13) 


In the rest of this section, we will present a method for solving (3.9)-(3.13). Consider a constrained 
optimization problem as 

max f(x), ifR", 
subject to gi{x) < 0 ,i = ... ,m 

then we define the i th constraint to be active (at a solution y) if gfiy) = 0. We will now show that in 
optimality, ( |3. 101 is active and the optimal solution always lies on the most restrictive constraint(s), i.e., the 
constraint(s) that impose the largest restriction on the dose that can be delivered to the tumor. 


(3.11)..(3.131 and furthermore constraint (3.10) is active in optimality. 


Lemma 1. The optimal X* and Y* in (3.9) lie on the feasible boundaries of the region defined by 


We provide the proof of this result in the appendix. Now we discuss the impact of adding X 2 < N„ 
and X 2 > Y on optimal solutions of (3.9). 


,Y 


Lemma 2. At the optimal solution to ( |3.9| ), the constraints X 2 < N max Y and X 2 > Y are either inactive 
or the optimal solution occurs at corners (X^, W 1 ^) and (X( Nrnax \Y( Nrnax 'l'). 

The proof of this result is provided in the appendix. As a direct result of Lemma [l] and Lemma [2] 
we know that the optimal pair (X*,Y*) lie along the feasible boundary of the convex hull CTLN mao; - For 
each feasible pair we can find the unique value z(X,Y) that gives equality in (3.10). Thus to solve the 
optimization problem we find the pair (X, Y) on the feasible boundary that maximizes the function i(X, Y). 


In order to solve problem (3.9) we can find Zi = f (Xj, Yf) for each 1 < i < fc, and then look at 


ma,x{zi — g(N) : 1 < i < k}. This process will tell us the optimal corner, but it does not necessarily tell 
us the optimal pair (X*,Y*). In order to find the optimal pair it is necessary to consider the edges of the 
polygon because it is possible to construct numerical examples where the optimal solution does not occur at 
a corner. Note that point (X i; Yj) in is only connected to the points (Xj_i, Yj_i) and (X i+ i, Y)+i), 

where (. X 0 ,Y 0 ) = (X k ,Y k ) and (X k+1 ,Y k+1 ) = (X^Y^.We therefore define the vector valued function for 
1 < i < k 

(X z (t),Y,(t)) = {tX i + (l-t)X i+1 ,tY i + {l-t)Y i+1 ), 0<t< 1, 
and the inverse function 


z i (t) = {z-g(\X i (t) 2 /Y i (t) 1) 




z-aX j(<) 
Vi(t) 


f(a,/3)d/3da = (1 -p z )P(a > 0,/3 > 0)}. 


(3.14) 


Note that the value of z in (3.14) for every pair (Xj(t), Y)(t)) can be computed using bisection method. 
The minimum number of treatment sessions for a given (X,Y ) is |"X 2 /T], thus g ([ Xfit) 2 /Yi{tf\) computes 
the reproduction effect for every (X, (t), YJf )). We use Algorithm [2] in appendix which is designed based on 
the branch and bound approach to find the optimal solution of (3.9). On each edge of feasible region, the 


branching is done on variable t and the global optimal solution is found via searching through all sub-optimal 
solutions on each edge. The choices for upper and lower bounds in each subproblem is given in the appendix. 
The optimal number of treatment sessions for an optimal pair of (X*,Y*) is N* = \{X*) 2 /Y*~\. 

The previous result shows how the solution (X*, Y*, N*) to the simplified versions of problems posed in 
Section 2 can be found. In E3, authors proved that (X, Y) transformation is indeed possible and is without 
loss of optimality. Moreover, we can derive optimal {d*,... ,d* N } with the following result. 
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Theorem 3. Optimal solution of {d\,... ,d* N } retrieved from (X*,Y*, N*) takes one of the following two 
forms. 


1. N*Y* = (X*) 1 2 3 for some N* £ {1,..., N max }. In this case optimal schedule is given by d* = X*/N* 

for i = ,N*. Note that if N* = 1 then the schedule is hypo-fractionated, and if N* > 1 then the 

schedule is hyper-fractionated. 

2. (X*) 2 /Y* is not an integer: In this case, the optimal solution given by the following. Choose a positive 
integer j less than (X*) 2 /Y* and set 


d* = -..=£ = 


jX* + \J(N* - j){jN*Y* — jX* 2 ) 


jN* 


d* — ■ 
a j+1 “ 


— d* — 

— a N , — 


X*-jdj 

N*-j 


(3.15) 


The proof of above result can be found in the appendix. 


4 Results 

In this section, we first discuss the effect of uncertainty on the structure of the optimal schedule. Then the 
application of nominal and robust optimization to the treatment of head and neck tumors via radiotherapy 
will be discussed. We will describe the data set and parameters that were used in our numerical experiments, 
then the solution to the nominal and robust optimum dosing schedules will be explored. At the end of this 
section the sensitivity of the optimal solution to model parameters is studied. 


4.1 Effect of uncertainty on the optimal solution 


Here we study the effects of parametric uncertainty by considering what happens to the optimal solution 
when the linear and stochastic robust formulations are used instead of the nominal formulations. We use 
(X n ,Y n ,N n ) and (X r ,Y r , N r ) to denote the optimal solutions to the nominal and robust (either stochastic 
or linear) problems. Throughout this subsection, we assume that the nominal and mean values of fli/oti and 
Si are equal to Pi/di and Si, respectively and the probability pi is greater than 50%. By imposing these two 
assumptions the feasible region of stochastic and robust problems becomes a subset of the feasible region of 
nominal problem (see Figure|2|. Note that the results presented in followings are valid only if we do not allow 
uncertainty in a and /3 (using same objective function as (2.3)). We will study the effects of uncertainty in 


a and /3 in the specific context of head and neck cancer in the next section. 


The feasibility region in the nominal formulation (2.3) is defined by several inequalities. Note that 


after using the transformation (3.7), these become linear inequalities. By introducing linear or stochastic 


uncertainty every line segment associated with each inequality will be broken down into two line segments 
with different slopes where each segment passes through ( Di , ^ ). There are different possibilities, depending 
on the amount of uncertainties in model parameters, slope of objective function (a and /3), reproduction rate 
of the tumor and tumor kick-off time, for how the optimal solution of the robust and stochastic problems 


relates to the optimal solution of the nominal problem. The feasible region of (2.51 and (2.6) can have either 
more, less or the same number of corner points as compared to the feasible region of (2.3) (see Figure |2]). 
We have three different scenarios for the robust or stochastic optimal schedule. 

1 . N r < N n : In this case we require that the total dose delivered to the tumor decrease, i.e., X r < X n . 
However we may have an increase or decrease Y r depending on the parameters (compare Figur^2]-A 
and[2}B). 

2. Ay > N n : Here we always have Y r < Y n . We can not say much about X r and it may be greater or 
smaller than X n (compare Figure [2]- A and[2]-C). 


3. N r = N n : In this case since the feasible regions of (2.5) and (2.6) are subsets of the feasible region of 


(2.3), we require that the total dose and total dose squared delivered to the tumor decrease or stay 


same. If for two normal tissues, i and j, we have Di = Dj, Ni = Nj and ( X n ,Y n ) = (Di, then 
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if the corner (X n ,Y n ) stays feasible after adjusting feasibility region based on model uncertainties, we 
will have X r = X n and Y r = Y n (see Figure [2]- A and[2]-D). 


To help illuminate the reasoning behind these ideas we give a short proof of (2). 

First note that the optimal solution satisfies the condition X 2 < NY. We know that the objective 
function is decreasing in N, thus the optimal N will be the smallest integer satisfying the constraint X 2 < 
NY, i.e. N n = |~(X*) 2 /Y*"|. 

Let to* = X 2 /Y* where (X*,Y*) are the optimal pair for the nominal problem and note that to* is 
not necessarily an integer. Assume N r > N n and then obviously N r > N n > to*. For any real number 
mi £ {N r — 1, N r \, let (Xi, Y\) be the intersection of line X 2 = miY with the feasible region of the nominal 
problem. Since TOi > N n and the slope of line segments constructing this feasible region are negative, then 
we have Y\ < Y*. Next define (X2W2) as the intersection of the line A' 2 = miY with the feasible region 
of the stochastic/robust problem. Since the feasible region of stochastic/robust problem is a subset of the 
feasible region of the nominal problem defined by line segments with negative slopes, then we have Y2 < Y\. 
Since it is true for every m\ € {N r — 1, N r \, then we can set toi = X 2 /Y r , and thus conclude Y r < Y\ < Y*. 
A similar argument shows that if N r < N n , then we have X r < X*. 

For model (2.5), the changes in X r and Y r depend on both N r and the amount of uncertainty in model 
parameters. Larger uncertainties (large 1/) result in larger reductions. In model (2.6), these reductions not 
only does depend on N r and the amount of uncertainty (defined by the variance of uncertain parameters) 
in the random variables, but also a key factor is the risk tolerance of the decision maker which is defined by 
Pi. If we have pi / 1, the reduction in (X r ,Y r ) in (2.6) is smaller than (2.5). 


4.1.1 Hyper vs Hypo-Fractionation 

One interesting question we can investigate is, when is a hyper-fractionated schedule preferable, and how does 
this compare to the setting without parameter uncertainty? In order to answer this question we make some 
simplifying assumptions. First we ignore tumor repopulation, and second we assume that there is tumor radio 
sensitivity parameters a and p are deterministic (i.e., known values). If pi > 0.5 and P{Si{Pi/ai) < 0) « 0, 
then we will have following different scenarios. 

• If {a/f. 3) < mini<i<M 1/6®, then a hypo fractionated schedule is optimal. 

• If {a/P) > maxix^jy/ 1/6', then a hyper fractionated schedule with N max fractions is optimal. 

• If mini<,<M 1/6® < {a/P) < maxi<i<M 1/6), either a hypo or a hyper or an unequal multiple dosage 
solution can be optimal. 

The constants bi and 6' are related to the distributions of /3j/cq and Si and are defined in the appendix. 

In the nominal setting in the absence of tumor reproduction, Ea showed the following. 

• If {a/P) < min-| < ? ;<m( d,/ Pi) /<5,;, then a hypo fractionated schedule is optimal. 

• If {a/P) > maxi<i<jv/(cL//3j)/<5j, then a hyper fractionated schedule with N max fractions is optimal. 

• If mmi<i<M(&i/Pi)/Si < {a/P) < maxi <i<Ai(ai/Pi)/8t, either a hypo or a hyper or an unequal 
multiple dosage solution can be optimal. 

Based on the previous two sets of results we see that we are using 1 /6, or 1/6' instead of {di/pi)/Si. 
Since pi > 0.5, it means 1 /6® < (a,//3,)/6j < 1/6', therefore we are making vam.i{{di/p/) / 5{\ smaller and 
ma,Xi{{di/Pi)/Si} greater. The result of this is that we increase the chance of having min,(l/6,) < {a/p) < 
max,; (1/6'), where either a hypo or a hyper or an unequal multiple dosage solution can be optimal. 


4.2 Application to head and neck tumors 

In order to estimate head and neck tumor radiobiologic parameters, we use the data set in [35] ■ To improve 
estimation accuracy, trials with the same properties such as total dose administered, number of fractions 
and treatment duration were merged. Model fit is carried out by minimizing the weighted error between the 
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model predictions of survival probability and the observed values of survival probabilities in trials associated 
with different schedules. In particular, the results of K trials have been considered. The trial outputs are 
survival fraction for each trial f\...., fx- We assume that the tumor cell population regrows exponentially 
after irradiation with a time lag of T k and rate 7 . Then following fractionated radiotherapy with X = 
and Y = y/( V 1 the population of tumor cells T units of time after start of therapy is given by 

N(T) = N{ 0) exp[-aX - jSY] exp[ 7 (T - T k )+], (4.16) 


As a simplification we say that recurrence is only detectable if N(T)/N( 0) > 1, i.e., if the tumor is bigger 
than its size at the start of therapy. There are K total trials and the total radiation in trial i is Xi and the 
sum of the doses squared in trial i is Y), and that the radiotherapy lasted for T ri days. Then from (4.16) we 
want to choose the distribution of a and /3 such that for each 1 < i < K 


P(exp[-aX - /3Y] exp[ 7 (T - T k ) + ] < 1) = fc 


Assume that the distributions of a and /3 are characterized by the joint density f(x',y,9) where 9 is a 
parameter that specifies the distribution. We define the probability that (a, 0) take values in some set as 
(with the 9 dependence explicit) 


Pe(a 1 < a < a 2 , h < Pi < b 2 ) 




f(x,y;9)dxdy. 


We assume that 9 takes values in the space 0. For each trial 1 < i < K we define function 


4>i{9) = P e (exp[-aX l - 0Yi] exp[ 7 (T - T k )\ < 1). 

Then our procedure for finding the best parameter set is to solve the minimization problem 

K 

min ^71^(60 - / 4 ) 2 

i= 1 

where rii is the number of patients in i th trial. Simulated annealing algorithm is utilized to find the optimal 
values of above model. We assume that the radiosensitivity parameters of LQ model, a and (3 are distributed 
based on two independent normal distributions with means y a and yp and standard deviations a a and erg, 
respectively. The reproduction rate 7 (= ^^) and Kick-off time, T k , for head and neck were selected to be 
0.003 per day and 21 days respectively [33]. The parameter T was set to be 5 years and the nominal values 
of a and /3 were set to be equal y a and yp. All parameters are summarized in Table [l] 

We consider six different normal tissues involved in the treatment of head and neck carcinomas ([?] and 
[33]). The nominal values and confidence intervals for p/a for various normal tissues were extracted from 
[35], [41], [10], [25], [30], [46], [24] and [31] and are listed in Table [ 2 J We assume that the ratio of {0/caY^i 1 
for normal tissues are distributed based on normal distributions with means y(p/ a ) i and standard deviations 
tj(p/ a )i- The values of means y(p/ a )i were set to the average of lower bound and upper bound of confidence 
intervals reported in above references. Also the standard deviation of {/3j/a,}^f^ associated with different 
normal tissues were computed based on their confidence intervals given in references. Normal distributions 
with parameters ySi and are considered for ■ Mandible and spinal cord are considered as serial 

structures and the data reported in [13] is utilized to compute their distribution parameters. Brain stem 
and parotid glands are assumed to be a serial and parallel tissues, respectively, and their parameters are 
estimated from data reported in [531 . In these papers, average and standard deviation of dose absorbed by 
a normal tissue for a given dose radiated to the tumor are reported. In 7], the values for planned dose, 
actually delivered dose and re-planned dose have been reported. We used these values to obtain a range for 
sparing factors for larynx and skin. Skin is considered a serial structure and larynx is considered a parallel 
structure. Nominal values of {Pi/aiY^zf 1 and sparing factors are set to the y((3/a)i and ys it respectively. 
The tolerance dose values for various normal tissues were computed from m, delivered in 35 fractions. We 
assume that patients may be treated at most in seven weeks and they visit the clinics three times a day, 


12 



n = 3. By considering 5 working days every week, we can compute the maximum number of allowable 
fractions as N max = 7x5 x n = 105. 

In order to understand the effects of parametric uncertainty on the structure of the optimal schedule, 
we consider one additional setting for /?. In particular we change the values of up and crp to 0.0001 (we 
call this scenario case 2 and the original values are case 1). The consideration of small values for /3 enables 
us to study the effect of uncertainty on schedules with a large number of fractions. The maximum dose 
constraints of 32 Gy for parotid glands results in optimal schedule with small values of total dose. Since in 
clinical practice larger dose has been used, we report optimal solution for two cases, including parotid glands 
in our constraint set and excluding it. Table [3] displays the optimum schedule for different models (optimal 
doses can be calculated from Theorem ([3])). Our numerical results for head and neck tumor show that 
the presence of uncertainty changes the optimal schedule to a schedule with larger dose delivered in more 
fractions (smaller total dose squared) for case 1 and almost same total dose delivered in smaller fractions 
for case 2 . As expected, the value of the objective function z* in (2.6) is decreasing in the probability of 
having the actual tumor BED less than the optimal value of z*. The z* value drops from 18.01 to 6.22 if we 
increase p z from 50% to 90%. 

Figure [i] plots the N* in (2.6) for different values of 7 = . The optimal value of the objective function 

is a non-increasing function in 7 . For short schedules, since the proliferation effect is negligible g(N) ss 0, 
the objective solution is robust to drifts in 7 . However in fast growing tumors, long treatment times have a 
negative effect on the treatment outcome. 


5 Conclusion 

In this work, we have analyzed the problem of finding optimal radiation administration schedules considering 
various types of normal tissues in the presence of model parameter uncertainty. In particular, we aimed to 
identify the optimized total dose, number of fractions, dose per fraction and treatment duration for a variety 
of formulations considering different types of uncertainty. We used the traditional linear quadratic model 
including tumor proliferation to investigate the dynamics of radiation response considering two uncertainty 
sets. First we assumed that only a range of possible values is known for the model parameters, ^ and sparing 
factors of normal tissue, <5, are known. We presented robust formulations of our optimization problem that 
are immune to realizations of the uncertain parameters so long as they lie within their respective ranges. 
Since using the ranges alone may lead to an excessively high level of conservativeness, in the second phase, 
we adjusted our formulations for the cases that uncertain parameters are distributed as continuous random 
variables with known probability density functions. Here we imposed the risk aversion factors in the objective 
function and the feasibility of constraints using some pre-dehned probabilities. 

We used the transformation introduced in Wl defining the total radiation as X and sum of doses squared 
as Y, and showed that our problem can be significantly simplified and easily solved in two dimensions when 
uncertainty in at and /3 are disregarded in the problem. In this case we observed that if the constraint 
X 2 = NY is active in optimality, then the largest possible BED to the tumor can be given in an equal- 
dosage schedule, otherwise the optimal solution is a semi-equal dosage schedule. When we consider a and (3 
as two continuous random variables, the problem becomes more challenging and optimal solution can happen 
at a non-corner point. In this case, first we have shown that the optimal value occurs at the boundaries of the 
feasible region defined by normal tissues BED constraints. We the designed a branch and bound algorithm 
to solve these stochastic models to optimality. Saberian et al. m have recently proposed a method to 
extract the optimal doses d\,d%, ... ,d^ given X*, Y*. However their approach fails if the optimal number 
of radiation sessions becomes larger than 2 and (X*) 2 /Y* is not an integer (note that d\ is not necessarily a 
positive real number in case 3 of Theorem 1 in (37] when {X*) 2 /Y* > 2). Here we showed that a semi-equal 
dosage schedule is optimal where d* = ■ ■ ■ = d* and d* +1 = ■ ■ ■ = d* N » for an integer j < x y . 

As a generalization of our results, we observed that when the presence of uncertainty does not change 
the structure of the optimal solution, it is preferred to administer same or smaller total dose and total dose 
squared. However if we have larger (smaller) treatment sessions in probabilistic or robust solution compared 
to nominal schedule, a reduction in total dose squared (total dose) will be seen. 

Using data gathered previously jM|, we parametrized the uncertainty in a and /3 to investigate the 
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behavior of optimal schedules for the head and neck tumors. For the numerical results, we assumed that the 
head and neck cancer site includes six normal tissues, spinal cord, brain stem, skin, mandible, larynx and 
parotid glands. The uncertainties in normal tissues have been estimated based on various data sets in the 
literature. The nominal optimal solution is a hypo-fractionated schedule changing to a schedule with larger 
total dose delivered in more fractions in the presence of parameter uncertainty. We found that when we 
consider small values of /3, the optimal schedule is a hyper-fractionated schedule with maximum allowable 
fractions. In this case the robust solution has an insignificant change in the optimal total dose and total 
dose squared for different schedules, however the optimal number of fractions decreases in some cases. We 
saw that as the tumor regrowth rate increases, shorter treatment are preferable. 

There are several possible extensions to this work that we plan to consider in the future. For example, 
this work does not incorporate spatial structure of the tumor, including possible spatial heterogeneities in 
the parameters a and j3. Another possible extension is the incorporation of repair effects, this would be 
useful if we wanted to consider shorter inter fraction periods. Lastly, it would be interesting to incorporate 
immune response and how inter-patient heterogeneity in immune response could impact the design of optimal 
fractionation schedules (see [ZD]). 
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6 Appendices 


6.1 Derivation of the robust reformulations 

6.1.1 Non-probabilistic robust formulation 


In this section, we derive a computationally tractable solution to the robust optimization problem (2.3). In 
the non-probabilistic robust formulation we do not allow any violation of the normal tissue constraints for 
any parameters taking values in the sets (2.4). Therefore the robust counter part of (2.3) associated with 
uncertainty sets defined in (2.4) is found by solving 


N 

max adn + /3dR — g(N) (6.17) 

d,->o,JVez+ ^ 3 0 ' 

3 =1 

subject to 


sup 






[— -h,— + k] and Si e [Si 
ai on 


hi , Si + hi] / < Di 


J 2 d j Vi. 


j =i 


Note that when the supremum happens when ^ and 6',; take their upper bounds in the sets 

(2.4), otherwise the supremum is achieved in lower bounds of — and Si defined in (2.4). We now replace the 


problem (6.17) by a formulation using the supremum of BED constraints: 


max 

dj>0,NeZ+ 


N 

ad J 

3 =1 


Pfi - g(N) 


(6.18) 


subject to 

where is defined as follows 


n B2 n • • • n Bm—i n Bm 


N 


N 


. i=1 3= 1 


Z^ d 3 


N 

CLj y d] < Dj 

i =1 


D 


N 








N, 


N 

E d > 

3 =1 


A 


i=i 


and a,i = (<5j + hji^r: + k) and a[ = (Si - hM^ - k)- 


6.1.2 Probabilistic optimization models 

We now describe a formulation that assumes that a, j3, and ^ and Si are random variables with known 
probability distributions, and we use our knowledge of this uncertainty to enforce the constraints in a 
probabilistic fashion. We require that the probability of violation of the BED constraint in i th OAR is at 
most 1 — pi. Written mathematically we have 


/ N 

p Y hh 


- 1 - A, 


(6.19) 


Note that from a biological point of view, it is impossible for a, /3, — and Si to take negative values, and 
since in literature it is often assumed that these parameters are normally distributed, we need to add the 
non-negatively conditions above. By knowing the cdf of — and Si (recall that we assume these random 
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variables are independent), we can easily compute the cdf of their product. Let Hi be the cdf of In 
order to satisfy |6 .19| for random variables 5i and we must have: 


N 


p ((Z! d2 o - ?r)ft— < A - ft-A— > 0 I > pA( o) v*. 

i=i A “* 3= 1 


A/ - 


.ft 


( 6 . 20 ) 


Note that we assume ft £ [0,1] and thus the event {A > 0, ft > 0} is equivalent to {ft^ 1 >0}. If we consider 

the events in (6.20) for the positive and negative values of YljLi dj — jf~, we can rewrite the condition (6.20) 
as where 


N 


Df 


N 


N 


s.={Y. d YtX d i+ b -T. d Y D - 


j=i 


i=i 


j=i 


A 

1 W 


JV 


D 


N 


N 


U 1 5Z ft + - A + 


. i =1 


i=i 


j=i 


A 

JVi 


and bi = H i X (1 - (1 -pi)Hi( 0)) and b\ = H i X (1 - PiHf 0)). 

We require that optimized schedules obtained by our robust formulation result in an objective 
value which exceeds level 2 , with a probability at least p z , i.e., 


N 


A =1 


a > 0,0 > 0 | >p z . ( 6 . 21 ) 

Note that for a fixed N, g(N ) is a constant and does not depend on di and therefore we can remove 
it from above probability. We can simplify (6.21) to get the constraint 


/*»/ 

Jo Jo 


’--<*T,y =1 dj 

A =1 dj f(a, (3)d(3da < (1 — p z )P(a >0,(3> 0). 


>o Jo 

6.2 Proof of technical Lemma [l] 


Proof. Assume z*,X* and Y* are optimal points to (3.9). If X* and Y* lie in the interior of 
the feasible region, then there exist AA' > 0 and AY > 0 such that the pair where 

X' = X* + A A and Y' = Y* + Ah', is a feasible solution. The left hand side of (3.10) is a 
decreasing function in X and Y. therefore we can replace (X*,Y*) with (X',Y') and increase z 
without violating feasibility of constraints set defined in (3.10) ... (3.13). Therefore there exists a 
feasible z which is strictly greater than z* and it contradicts the assumption that z* is an optimal 
solution to our problem. Therefore the optima must lie on the boundaries of the feasible region. 
Also if (3.10) is not active in optimality, we have 


/*/ 


z — cxX 
Y 


f(a, f3)d/3da < (1 —p z )P(a > 0,j3 > 0) 


and we can increase z without leaving feasible region which contradicts the optimality assumption. 

□ 


6.3 Proof of technical Lemma [2] 


Proof. Let CT-L^max be the feasible region of (3.9) defined by (3.11), (3.12) and (3.13). Based on 
lemma [I] the optimal solution of (3.9) lies on the feasible boundaries of C 7 v maa ,. If X 2 <N max Y and 


16 






















X 2 > Y are redundant constraints ((X*,Y*) obtained by ignoring these constraints satisfy these 
constraints), then X 2 < N max Y and X 2 > Y are inactive constraints in optimality. Otherwise, 
consider the three corners pi = (0,0), p 2 = (X^ 1 ), Y^) and p 3 = (X < ' Nrnax \Y l ' Nrnax ')). As we move 
from pi toward P 2 or from p\ toward P 3 , we can increase both X and Y. At the end of two line 
segments pop{ and P 0 P 2 , we are at a feasible solution (p 2 or p 3 ) with maximum X and Y, and since 
the objective function in (3.9) is increasing in both X and Y, we see that the maximal value of 2 : 
is obtained at either p\ or P 2 . □ 


6.4 Proof of Theorem [3] 


Proof. To establish the result in case 1, we can easily check that equation (3.7) holds for 
d\,..., djy* = X*/N* if N*Y* = (A'*) 2 . In the second scenario, from straightforward calcula¬ 
tions, we observe that there is always a solution to our problem in the following form: 


d * 1 = ■■■ = d* = d, d * j+1 = ■ ■ ■ = d* N , = w. 


We can now solve for w and d in (3.7), and establish that 


w = 


X*-jd 
N* — j 


and 


_ jX* + i/j' 2 (X *) 2 - jN*((X *) 2 - (N* - j)Y*) _ jX* + VCiV* - j)(jN*Y* - j(X*) 2 ) 

jN* ~ jN* 

It then remains to establish that d and w are non-negative real numbers. First observe that we 
require that j < {X*) 2 /Y* < N*. It follows from this that d is a positive real number, and thus 
w is a real number as well. It then remains to establish that w is non-negative. This is of course 
equivalent to showing that X* — jd > 0. Note that 


X*-jd 


(N* - j)X* - y/j(N* -j)(N*Y* -(X*) 2 ) 

(N* - j)X* - yJjY*{N* — j ) (N* - (X*) 2 /Y*) 


and therefore 

( Y*l 2 

X*-jd>0^X ¥ L>j. 

The result then follows from our conditions on the integer j. □ 
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6.5 Algorithms 

6.5.1 A Feasible Region Creator Algorithm 


Algorithm 1 Constructing the corners of feasible region defined by (3.8), X 1 > Y and X 1 < NY 


D- D■ D- 

1: Define C\ as the set of all line segment connecting (0, + ~ff:) an d (Dj, jfr) and £2 as the set of all line segment passing 

through (Di , and (A + cl 2*-, 0) for i = 1,... , M. 

2: Use Bentley-Ottmann algorithm [3] for listing all crossings in the set of {£1 U £ 2 }, call it £. 

3: Let £ be £ U {(0, + ^r~)} U {(A 2 + c i 2 jv“- 0 )} where i\ = argmin k {^ and i 2 = argmin k {D k + c' k ^-} 


N k - 


k N k - 


4: Compute V = {(??i, ?;])...., (v p , v ' v )}. as the set of all pairs in £ satisfying (|3.S[l for all i = I..... M, sorted in increasing 


order by the x-coordinate. 

5: CH\ {(0, 0)} 

6: For i <— 1 to p — 1 

7: Compute ( Xi,yi ) as the intersection of line passing through (?; v ,?/) and (??,—i,??] ]) arid y = x 2 . 

8: If Vi < Xi < Ui-|-i and yi > 0 

9: index[l]<— i 

10 : CHi^CH-i U{(xi, yi )} 

11: Break 

12: End If 

13: End For 

14: For i «— 2 to Al maa; 

15: Off,; C'A-i 

16: For j index]* — 1] to p 

17: If v 2 < iv'j 

18: 

19: indexfil <— j 

20: Else 

21: if j ==index[i — 1] then index)?] <— index]* — 1], else index]?] <— j — 1 

22: Break 

23: End If 

24: End For 

25: End For 

26: For * <— 2 to N, nax 

27: Let (»<,{/,) be the intersection of x 2 = **/ and line passing through Oi ndex [i] > *4dex[i]) and 6Wex[i]+l> u i ndex [i]+l) 

28: CF^OT.Ufxi,!,,) 

29: End For 


6.5.2 Branch and bound algorithm 


Note that for every line segment defined by two points (X \, Yi) and (X 2 , Y 2 ), we can compute 
the lower bound using 

®lb = max{zi(0),zi(l)} 

and compute the upper bound using 


*&ub — 


{z-g( r 


min(Ni, X 2 ) 
max(Yi, Y 2 ) 


1 



z — a max(Y ^ , JNT 2 ) 
max(Y 1 ,Y 2 ) 


f(a, /3)d/3da = (1 — p z )P(a > 0, (3 > 0)}. 


On each edge, the branching variable is t € [0,1]. At node i, we store the lower bound and upper 
bound of t as lbt(i) and ubt{i). Similarly lower and upper bounds of optimal solution at node i are 
stored in lb z (i ) and ub z (i). Note that a(i) = {0,1} indicates the state of the node in our optimization 
tree, if a(i) = 1 then our node is considered as active node and further partitioning can be proceeded 
through that branch, otherwise we consider that node as an inactive node. Branching in node i 
continues in this manner until there are no active nodes in that branch or ub z (i) — ub z (i ) < e, where 
e is the given accuracy for optimal value. Since the objective value of an optimal solution cannot be 
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Algorithm 2 Branch and Bound algorithm for maximization of (3.9) 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 


For each (Xi,Yi) E ^^Nmax » compute Zi( 1) and set z* = Zj( 1) and (X*,y*) = ( Xj,Yj ) where j = argmax^ z% 
For l = 1 to |C'H_/v rnaa . | — 1 Do: 
i <— 1. 

lbt(i) 0, ubt(i) «— 1. 

lb z (i) = max{ 2 ?(Xj( 0 ),yi( 0 )),«(X|(l),yi(l))}, ub z (i) = z(max{X z (0), X z (l)}, max{y z (0), y z (l)}). 
a(i) <— 1. 

While jy a(i) > 0 Do: 
ind=find(z|a(z) > 0). 

For j = l:length(ind) Do: 

i «— i + 1. 

lb t (i) = (/6 t (ind(j)) + ubt(md(j)))/2, ub t (i ) = u6t(ind(j)). 

Compute 


lb z (i) = max{ 2 :(X z (Z 6 t(i)),y z (/ 6 t(i))), 2 :(X z (w 6 t(i)),y z (?x 6 t(i)))}, ub z (i) = z(max{Xi(lb t (i)), Xi(ub t (i))}, max{y z (Z 6 £ (i)), Yi(ubt(i))}). 


13: If ( ub z (i ) — lb z (i) > e), then a(i) «— 1, else a(i) «— 0. 

14: /6t(ind(j)) = lb t (ind(j)), ub t (ind(j)) = (lb t (md(j)) + ubt(md(j)))/2. 

15: Repeat steps 13 and 14 with ind(j) instead of i. 

16: Update z* and (X*,Y*) for any lb z (j) > z*. 

17: End For Loop. 

18: For every j such that ub z (j) < z*, a(j) <— 0. 

19: Sort [lbt(:), ubt(:), lb z (:), ub z (:), a(:)] T column-wise based on lbt(i). 

20: End While Loop. 

21: End For Loop. 


smaller than a lower bound, active nodes with upper bounds smaller than an existing lower bound 
can be safely deleted (step 19). 

Remark 1. Algorithm [S] converges and terminates with certificate proving e-suboptimality. 
Number of line segments in partition is k. Note that total length of these line segments is 

-h(hlinitial); SC) 

r/^\ ^ -h(£2i n ;tial) 

nun L Q < --- 

Qefife k 

and hence for big k, at least one line segment has small length and having small length will imply 
that ub z (k ) — lb z {k) is small. 
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6.6 Figures 


Y Y 



Figure 1: Feasible region of optimization problems posed in section 2. A) This plot shows the 

feasible region of model (2.3). For each N, the optimal solution occurs at one of the corners of 

feasible region. B) This plot shows the effect of uncertainty in model parameters on feasible region. 

Every line segment is broken down into two line segments with different slopes and passing through 
£)2 

(. Di , -rfr). If a and /3 are fixed, the optimal solution for every fixed N lie on one of the corners of 
CHn■ Having a and /3 as random variables, the pair (X*,F*) can be located by searching on all 
line segments connecting (A' 1 ,!' 1 ) and (X Nmax , Y Nrnax ). 


20 












Figure 2: This plot shows the effect of uncertainty on the optimal solutions in the presence of 
three OARs. In each subfigure the feasible region is shaded. In panel (A) we plot the feasible 


region and optimal solution in the nominal case. Note we also label the points ( D[ 


D' 3 , 


( D ' a ) 2 


( D \) 2 

Ni 


and 


as alternative maximum tolerable doses. In panel (B) we consider a scenario where 

In panel (C) 


3’ T/vT 

the solution to the stochastic problem results in fewer total fractions, i.e., Nb < N n 
we consider different distributions for our uncertain parameters and we have a scenario where the 
optimal number of fractions in the stochastic problem is greater than the number in the nominal, 
i.e., Nc > N n . Finally in panel (D) we use the alternative maximum tolerable doses (d[, 


and and construct a scenario where the optimal number of doses is unchanged by 

parameter uncertainty, i.e., Nb = N n . If the optimal number of radiation sessions stays the 
same in the presence of uncertainty, it is required to deliver equal or less doses in the presence 
of uncertainty ( D ). If we have fewer (more) treatment sessions in probabilistic or robust solution 
compared to nominal schedule, a reduction in tSllal dose (total dose squared) will be seen (B and 
C). 











Sensitivity analysis of optimal number of fractions with respect to y 



Figure 3: This plot shows the sensitivity of treatment session in (2.6) assuming Pi = p z = 95% with 
respect to tumor growth rate 7 assuming T\ = 7 days. For short schedules the objective solution 
is robust to drifts in 7 and for fast growing tumors, long treatment time have a negative effect on 
the treatment outcome. 
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6.7 


Tables 


Parameters 

Values 

unit 

Mo; 

0.1708 

1 /Gy 

M/3 

0.0537 

1 /Gy 2 


0.2142 

1 /Gy 

vp 

0.0812 

1/Gy 2 


Table 1: Head and neck tumor parameters used for finding optimal schedule 


Parameters 

Spinal Cord 

Brain Stem 

Skin 

Unit 

M(/3/a) 

0.48 

0.39 

0.12 

1/Gy 

a (P/a) 

0.09 

0.05 

0.01 

1/Gy 

(/3/a - /,/3/a + /) 

(0.30,0.67) 

(0.30,0.48) 

(0.09,0.14) 

1/Gy 

PS 

58.52% 

74.92% 

25.29% 


os 

2.78% 

3.45% 

0.75% 


Di 

47 Gy 

50 Gy 

55 Gy 

Gy 

Parameters 

Mandible 

Larynx 

Parotid glands 

unit 

M(/3/a) 

0.46 

0.66 

0.24 

1/Gy 

a (P/a) 

0.05 

0.30 

0.07 

1/Gy 

(/3/a - /,/3/a + /) 

(0.36,0.56) 

(0.07,1.25) 

(0.10,0.38) 

1/Gy 

PS 

67.59% 

99.12% 

40.45% 


os 

10.56% 

0.07% 

0.97% 


Di 

60 Gy 

70 Gy 

32 Gy 

Gy 


Table 2: Normal tissue parameters 




With parotid glands 

Without parotid glands 

Parameters 

Formulation 

N* 

X* 

Y* 

Tumor BED 

N* 

X* 

Y * 

Tumor BED 


Nominal 

1 

13.5 

182.36 

12.10 

1 

13.5 

182.36 

12.10 

Case 1 

Robust 

2 

14.26 

139.51 

9.93 

35 

47.00 

63.11 

11.42 


Stochastic 

3 

18.34 

124.23 

4.36 

11 

31.86 

95.94 

5.08 


Nominal 

105 

33.79 

10.87 

5.69 

105 

56.26 

30.15 

9.53 

Case 2 

Robust 

62 

32.47 

17.01 

5.53 

105 

52.82 

26.57 

8.95 


Stochastic 

54 

32.48 

19.80 

0.90 

92 

52.84 

30.59 

1.41 


Table 3: Optimal solution to problems (2.3), (2.5) and (2.6) assuming pi = p z = 95%. In case 
1 we assume p a = 0.1708, pp = 0.0537, a a = 0.2142 and erg = 0.0812 and for case 2 we have 
p a = 0.1708, pp = 0.0001, a a = 0.2142 and erg = 0.0001. 
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